A Distance Measurement to M33 Using Optical Photometry of Mira Variables

Mira variable stars (hereafter Miras) are red giants located on the asymptotic giant branch (AGB) at their late stage of evolution. Miras are long-period variable stars, with periods ranging from hundreds to thousands of days, and they have large amplitude variations in the optical and near-infrared (NIR) bands (for examples, see Soszyński et al. 2009; Riebel et al. 2010; Whitelock 2012). Miras can be categorized as oxygen-rich (O-rich) or carbon-rich (C-rich) stars according to the nature of molecules that dominate their spectra; this depends on the C/O ratio presents in their atmosphere (Merrill 1960; Cioni et al. 2001; Riebel et al. 2010). Alternatively, division of Miras into O-rich or C-rich candidates can also be done by using photometric data when the spectroscopic data is lacking. For examples, Miras were categorized on the basis of period and Wesenheit indices W_I in Soszyński et al. (2005). On the other hand, Soszyński et al. (2009) found that both categories of stars can be separated using a (V − I) versus (J − K_s ) and W_JKs versus W_I diagram. In Yuan et al. (2017a) and Yuan et al. (2018), the authors also divided stars into subclasses by using a (J − K_s ) versus (H − K_s ) approach. In Lebzelter et al. (2018), the authors used the $({W}_{\mathrm{RP}}-{W}_{{K}_{s}})$ versus K_s to divide the O-rich and C-rich AGB stars based on the Gaia and 2MASS data. O-rich and C-rich Miras were also found to belong to two distributions in a period and (J − K_s ) diagram (Iwanek et al. 2021).

Since the first period–luminosity (PL; also known as the Leavitt Law) relation for Miras was published in Glass & Evans (1981) using NIR data, a number of researchers have attempted to derive and calibrate the Mira P–L relations (e.g., see Feast 1984; Feast et al. 1989; Kanbur et al. 1997; Whitelock et al. 2008; Yuan et al. 2017b; Bhardwaj et al. 2019; Iwanek et al. 2021; Ou & Ngeow 2022). The dispersion of the Mira P–L relation at maximum light was found to be smaller than that of their counterparts at mean light (Kanbur et al. 1997; Bhardwaj et al. 2019; Ou & Ngeow 2022).

Recently, we have derived and calibrated the I-band PL relations for Miras located in the Magellanic Clouds (Ou & Ngeow 2022). Hence, the main goal of this work is to test the applicability of our PL relations in measuring distances to nearby galaxies. We selected M33 for our test because there is a sizable sample of Miras found in M33 (Yuan et al. 2018), M33 is close enough such that light curves for these Miras can be retrieved from archives (Section 2), and there are numerous independent distance measurements to M33 so we can compare our derived distances with these independent measurements. The sample of 1781 M33 Miras compiled in Yuan et al. (2018) was classified into 88 C-rich Miras, 1265 O-rich Miras, and 428 of unknown type. Multiple reasons were accounted for the unknown type, as described in Yuan et al. (2018) and will not be repeated here. Using independent optical-band light curves we collected in Section 2, we redetermined the periods for this sample of M33 Miras, as well as their magnitudes at mean and maximum light, and reclassified the Miras with unknown type using a machine-learning approach. The results of this analysis are presented in Section 3. Since C-rich Miras could exhibit a significant long-term variation in their light curves (Iwanek et al. 2021; Ou & Ngeow 2022), which will affect the result of derived PL relation, we only use O-rich Miras (including those being reclassified) to determine the distance to M33 using our derived PL relation (Ou & Ngeow 2022) and compared to other distance measurements in Section 4, followed by the conclusion of this work in Section 5.

We collected the long-term (2003–2021) optical gri-band light curves for the 1781 M33 Miras compiled in Yuan et al. (2018) from various sources. We cross-matched this sample of Miras to the variable sources detected in the M33 variability survey (Hartman et al. 2006; hereafter HBS 2006), a time-series gri-band survey carried with the MegaCam mounted on the Canada–France–Hawaii Telescope (CFHT) from 2003 to 2005 (for 27 nights). Sparse gri-band light curves data were also collected (whenever available) from the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS; Kaiser et al. 2010; Chambers et al. 2016) Data Release 2. We then extracted the R_PTF-band light curves for the M33 Miras from the Palomar Transient Factory (PTF; Rau et al. 2009; Law et al. 2009) and the intermediate PTF (iPTF; Kulkarni 2013). These R_PTF-band light curves were calibrated to the r-band using the Pan-STARRS photometric catalog. Specifically, we performed differential photometry by selecting a number of suitable reference stars around each Miras, where the r-band magnitude for these reference stars are available from the Pan-STARRS photometric catalog. Together, the Pan-STARRS and the PTF/iPTF light curves spanned from 2009 to 2017. Finally, the gri-band light curves data after 2017–2021 were collected from the Zwicky Transient Facility (ZTF; Bellm et al. 2019; Masci et al. 2019; Graham et al. 2019; Dekany et al. 2020) Data Release 10 and the ZTF collaboration survey data. ⁶ Altogether, we collected gri-band light curves data for 1367 M33 Miras (see Table 1). An example of the collected light curve is shown in Figure 1. We noted that since the telescopes used in Pan-STARRS and PTF/iPTF/ZTF have an aperture of 1.8 m and 1.2 m, respectively, limiting magnitudes from these surveys are around ∼19.5 mag to ∼21.5 mag. In contrast the deep CFHT observations for the M33 variability survey can reach to a depth of ∼25 mag in r-band.

Zoom In Zoom Out Reset image size

3.1. Periods Determination

Since Miras exhibit a large amplitude variation in the optical bands, it is possible that for Miras located in a distant galaxy, only a portion of the optical-band light curve (i.e., around the maximum light) brighter than the limiting magnitude of a given survey can be detected. This is indeed seen in the light curves for M33 Miras collected from Pan-STARRS, PTF/iPTF, and ZTF, as demonstrated in Figure 1. On the other band, the deep HBS 2006 observations can sample the full-amplitude light curves, including the portion of the light curves around the minimum light. These combinations provide us an opportunity to test the period determination for distant Miras when only a portion of the light curves above a given detection limit is available.

We first determined the periods for our sample of M33 Miras using the multiband Lomb–Scargle (LS) periodogram, described in VanderPlas & Ivezić (2015) and implemented in the gatspy package, on the full set of gri-band light curves (whenever available). Errors on the determined periods were estimated based on bootstrap resampling method. We then checked the LS periods using a multiband phase-dispersion minimization (PDM) periodogram developed in Lee et al. (2021). If both periods agree (e.g., the periods are within 10% of each others) we adopted the LS periods, else we visually inspected the phased light curves and selected the period that resulted a smoother light curve. Upper panel of Figure 2 presents an example of the LS periodogram for HBS 2006–40671 (see Figure 1 for the observed light curves), at which a period of 652.7 ± 23.3 days was identified. ⁷ The periodogram from the multiband PDM, as presented in the middle panel of Figure 2, also picked up the same period as the LS periodogram. Even though the LS periodograms could have multiple peaks with similar heights, we applied both multiband LS and PDM methods for cross-check and validations to ensure the most probable periods were selected.

Zoom In Zoom Out Reset image size

Our determined periods are given in Table 2. In general, our determined periods agreed with periods presented in Yuan et al. (2018), as evident in the upper panel of Figure 3, validating our period-search method. For this sample of M33 Miras, only ∼4% of the Miras display a significant difference in the periods we found here and in Yuan et al. (2018), an example is presented in the bottom panel of Figure 3. Figure 4 shows that for these ∼4% of Mira, the overall phase dispersions for light curves folded with our determined periods are generally smaller than those using the periods from Yuan et al. (2018).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. Observed Properties of Miras in M33

ID	P (days)	σP (days)	Type	〈i〉	σi	iMAX	${\sigma }_{i\max }$	〈r〉	σr	rMAX	${\sigma }_{r\max }$	〈g〉	σg	gMAX	${\sigma }_{g\max }$	E(B − V)
01321450+3019349	260.7	3.42	O	20.82	0.06	20.44	0.06	22.23	0.07	20.83	0.08	24.48	0.19	23.95	0.07	0.054
01321654+3025260	304.1	7.72	O	21.24	0.04	20.64	0.12	24.19	0.03	22.95	0.2	25.65	0.02	24.18	0.21	0.051
01321897+3031226	255.4	8.24	O	21.7	0.11	20.53	0.81	23.91	0.06	21.66	1.89	25.69	0.13	23.89	0.03	0.05
01322179+3034063	356.0	0.24	O	21.05	0.1	20.4	0.36	23.69	0.04	22.75	0.99	25.28	0.13	24.37	0.37	0.049
01322351+3030590	252.3	0.13	O	21.58	0.14	21.01	0.5	24.52	0.09	22.77	1.03	27.05	0.15	24.49	0.32	0.049
⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯

Note. The entire Table is published in its entirety in the machine-readable format. A portion is shown here for guidance regarding its form and content.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

The deep but short duration (probably cover ∼1 to ∼2 pulsation cycles of Miras) of CFHT observations from HBS 2006 and the shallow but longer duration (covering more than two pulsation cycles) from the Pan-STARRS/PTF/iPTF/ZTF surveys represent two typical cases encountered in dedicated deep time-series observations (such as targeting a particular galaxy) or synoptic sky surveys. We tested two scenarios on period search by excluding the deep HBS 2006 light curves (and only using light curves from Pan-STARRS/PTF/iPTF/ZTF) and only using the HBS 2006 light curves. Comparisons of the periods found in these two scenarios to the assumed true period P using the full HBS 2006 and Pan-STARRS/PTF/iPTF/ZTF light curves are presented in Figure 5. Our test results show that to recover the periods, it is important to sample the full-amplitude light curves for Miras rather than only sampling a portion of the light curves around maximum light covering few pulsation cycles.

Zoom In Zoom Out Reset image size

3.2. Magnitudes at Mean and Maximum Light

3.3. Reclassification of M33 Miras with Unknown Type

In accordance with the results presented in Iwanek et al. (2021), we examined the period–color (PC) relations for the M33 Miras in various colors. Figure 6 reveals that the O-rich and C-rich Miras exhibit different distributions in the (J–K_s ) and (H–K_s ) PC relations, but their distributions were similar in the (i–H) PC relations. This implies that the colors of O-rich Miras is markedly different from that of C-rich Miras in NIR but not in the optical. That is caused by the C-rich Miras having higher abundance of circumstellar dust, which can significantly influence the near-infrared radiation (Iwanek et al. 2021; Ou & Ngeow 2022). Figure 6 also shows that the O-rich and C-rich M33 Miras can be well separated in the (J–K_s ) PC relation.

Zoom In Zoom Out Reset image size

The sample of Miras in the Large Magellanic Cloud (LMC; as compiled in Ou & Ngeow 2022), the O-rich and C-rich Miras have different distributions in the (J–K_s ) PC relation, as indicated in Figure 7. This suggested the O-rich and C-rich Miras can be classified in the (J–K_s ) PC plane via machine-learning (ML) techniques. We tested four ML classifiers, namely the perceptron learning algorithm (PLA), the logistic regression algorithm (LRA), the K-nearest neighbor (KNN) algorithm, and the support vector machine (SVM) algorithm.

Zoom In Zoom Out Reset image size

In the PLA, the training data were applied to identify a linear function that can categorize data into two groups; this linear function was used as the decision boundary. The LRA is similar to the PLA, but it uses a sigmoid function as an activation function. The KNN algorithm was applied to monitor the position of data in the feature plane and then investigate the types of nearest k neighbors. A final decision on Mira type was made based on the type with the largest number of neighbors. In this study, we set k = 15. The SVM algorithm was applied to identify a separating hyperplane for use as the decision boundary.

In our study, 50% and 50% of the Miras in the sample were used as training and test data, respectively. We have tested these machine-learning classifiers separately on the LMC and M33 samples. In both samples the accuracy were greater than 95%, however the decision boundaries were different. Hence, we adopted the M33 sample as the training data. The decision region obtained using each algorithms is presented in the left panels of Figure 8, and the confusion matrix derived from the test data is displayed in the right panels of Figure 8. All of the true–true ratios obtained from the test data were greater than 95%, except for those derived from the PLA; therefore, the relationship between the periods and (J − K_s ) colors could be applied to classify O-rich and C-rich Miras. We adopted the KNN algorithm to classify the 344 unclassified Miras presented in the Table 2 of Yuan et al. (2018). The number of O-rich and C-rich Miras was found to be 310 and 34, respectively; the corresponding results are presented in Figure 9 and Table 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Using the periods and magnitudes derived in Sections 3.1 and 3.2, the distance to M33 can be determined by fitting the PL relation given in Ou & Ngeow (2022). Before fitting the PL relation, we verified that the dispersion of PL relation is smaller at maximum light for the Miras in M33. We used a quadratic model to fit the P–L relation to the O-rich Mira i-band magnitudes at both maximum light and mean light (after corrected for extinction), where the quadratic model is given as:

$$\begin{eqnarray}&&m=a+{b}_{1}{(\mathrm{log}P-\mathrm{log}{P}_{b})+{b}_{2}(\mathrm{log}P-\mathrm{log}{P}_{b})}^{2},\end{eqnarray} \tag{ 1 }$$

with the break period adopted at P_b = 300 days (Bhardwaj et al. 2019; Ou & Ngeow 2022). The relevant results are presented in Figure 10. The PL dispersion was found to be 0.42 at maximum light and 0.57 at mean light, confirming the earlier results.

Zoom In Zoom Out Reset image size

When fitting the P–L relation to derive the distance to M33, we only used the quadratic model of the P–L relation, based on the LMC O-rich Miras, taken from Table 2 of Ou & Ngeow (2022), and fit to the O-rich Miras in M33 (including those reclassified in Section 3.3). A LMC distance modulus of μ_LMC = 18.49 mag (Pietrzyński et al. 2019) was adopted when measuring the distance modulus of M33. Since the P–L relation presented in Ou & Ngeow (2022) is in the I-band, we transformed our i-band magnitudes, together with the corresponding gr-band magnitudes, at maximum light to I-band using the transformation given in Tonry et al. (2012). We did not transform the i-band magnitudes at mean lights because the gr-band mean magnitudes are not reliable. The fitted P–L relation is presented in Figure 11, and we derived ${\mu }_{\max }^{M33}=24.67\pm 0.06$ mag at maximum light. Error on ${\mu }_{\max }^{M33}$ is the quadrature sum of the errors on μ_LMC (0.048 mag), the estimated error in the fitted P–L relation (0.022 mag), and the dispersion from the photometric transformation (0.017 mag).

Zoom In Zoom Out Reset image size

Our derived distance modulus falls in the middle of the range of the distance modulus found in the literature (see Table 3) and is the same as the recommended value of 24.67 ± 0.07 mag from de Grijs & Bono (2014). Yuan et al. (2018) found a distance modulus of 24.80 mag, albeit using a similar sample of O-rich Miras. We emphasize that both work used totally different data sets (such as observed in different filters) and methodologies (such as fitting of the P–L relations at mean and maximum light). If we adopted the periods derived in Yuan et al. (2018) and repeat the analysis, we found a distance modulus of ${\mu }_{\max }^{M33}=24.69\pm 0.06$ mag, which is fully consistent with the distance modulus if using the periods we found in Section 3.1.

Table 3. Distance Modulus for M33 from the Literature

μ	Method	Sample	Filter	Literature
24.50 ± 0.06	TRGB		g,i	McConnachie et al. (2004)
24.52 ± 0.19	Leavitt law	Cepheid	V,I	Lee et al. (2002)
24.53 ± 0.11	Leavitt law	Cepheid	B,V,Ic	Scowcroft et al. (2009)
24.55 ± 0.28	Leavitt law	Cepheid	V,I	An et al. (2007)
24.57 ± 0.06	JAGB		J,K	Zgirski et al. (2021)
24.57 ± 0.05	TRGB		i',g'	Conn et al. (2012)
24.58 ± 0.10	Leavitt law	Cepheid	V,I	Freedman et al. (2001)
24.60 ± 0.30	TRGB		B,V,R,I	Wilson et al. (1990)
24.62 ± 0.06	Leavitt law	Cepheid	J,H,K	Bhardwaj et al. (2016)
24.62 ± 0.07	Leavitt law	Cepheid	J,K,V,I	Gieren et al. (2013)
24.64 ± 0.09	Leavitt law	Cepheid	B,V,R,I	Freedman et al. (1991)
24.64 ± 0.15	TRGB		V,I	Galleti et al. (2004)
24.67 ± 0.05	JAGB		J,H,K	Lee et al. (2022)
24.67 ± 0.06	Leavitt law	Mira	I	This work
24.69 ± 0.07	TRGB		V,I	Tiede et al. (2004)
24.70 ± 0.11	TRGB		I	de Grijs & Bono (2014)
24.71 ± 0.08	JAGB		J,H,K	Lee et al. (2022)
24.71 ± 0.04	TRGB		F555W,F606W,F814W	Rizzi et al. (2007)
24.71 ± 0.04	Leavitt law	Cepheid	J,H,Ks	Lee et al. (2022)
24.72 ± 0.14	TRGB		V,I	Brooks et al. (2004)
24.72 ± 0.11	TRGB		J,H,K,I	Lee et al. (2022)
24.76 ± 0.02	Leavitt law	Cepheid	B,V,I	Pellerin & Macri (2011)
24.80 ± 0.06	Leavitt law	Mira	J,H,Ks	Yuan et al. (2018)
24.80 ± 0.10	Leavitt law	Cepheid	B,V,R	Metcalfe & Shanks (1991)
24.81 ± 0.04	TRGB		V,I	Lee et al. (1993)
${24.81}_{-0.07}^{+0.19}$	TRGB		F555W,F814W	Kim et al. (2002)
24.82	TRGB		V,I	Salaris & Cassisi (1997)
24.84 ± 0.10	TRGB		F814W,F606W	Vivian et al. (2009)
24.93 ± 0.18	TRGB		B,V,I	Ferrarese et al. (2000)

Download table as:  ASCIITypeset image

In this work, we aimed to derive the distance to M33 using Mira variables. We started by redetermining the pulsating periods for 1378 Miras in M33 using available gri-band light curves. While performing the period analysis, we found that it is crucial to sample the full-amplitude light curve (in optical bands) rather than portions of the light curve around maximum light. This will be particularly important for Miras discovered in distant galaxies by the Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST; Ivezić et al. 2019), or other similar sky surveys, where it may not be possible to sample the light curves around minimum light due to limiting magnitudes.

In addition to period analysis, we showed that O-rich and C-rich Miras can be separated on (J − K_s ) PC plane using ML techniques, and hence we reclassified those Miras with unknown types in Yuan et al. (2018). Using all available O-rich M33 Miras, we demonstrated the P–L relation has a smaller dispersion at maximum light. Finally, we derived the distance modulus to M33 after transforming the photometry to I-band and fitted with the quadratic P–L relations given in Ou & Ngeow (2022). The derived distance modulus is 24.67 ± 0.06 mag using the P–L relations at maximum light, which is in good agreement with literature values. Our work demonstrated that both P–L relations at maximum and mean light for Miras can be used in distance scale measurements, this will be very useful in the era of LSST.

We thank an anonymous referee for the suggestions to improve the manuscript. We are thankful for funding from the Ministry of Science and Technology (MoST, Taiwan) under the contract 107-2119-M-008-014-MY2, 107-2119-M-008-012, 108-2628-M-007-005-RSP, and 109-2112-M-008-014-MY3. A.B. acknowledges funding from the European Unions Horizon 2020 research and innovation program under the Marie Skodowska-Curie grant agreement No. 886298.

Based on observations obtained with the 48-inch Samuel Oschin Telescope at the Palomar Observatory as part of the Zwicky Transient Facility project. ZTF is supported by the National Science Foundation under grants No. AST-1440341 and AST-2034437 and a collaboration including current partners Caltech, IPAC, the Weizmann Institute of Science, the Oskar Klein Center at Stockholm University, the University of Maryland, Deutsches Elektronen-Synchrotron and Humboldt University, the TANGO Consortium of Taiwan, the University of Wisconsin at Milwaukee, Trinity College Dublin, Lawrence Livermore National Laboratories, IN2P3, University of Warwick, Ruhr University Bochum, Northwestern University and former partners the University of Washington, Los Alamos National Laboratories, and Lawrence Berkeley National Laboratories. Operations are conducted by COO, IPAC, and UW.

The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen's University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation.

Software: astropy (Astropy Collaboration et al. 2013, 2018), gatspy (VanderPlas & Ivezić 2015).